Vamos a ver que existen varios escenarios de raÃces unitarias en modelos ARIMA
library(urca)
library(forecast)
library(tseries)
library(lmtest)
library(uroot)
library(fUnitRoots)
library(aTSA)
######Ejercicios de Simulación#####
###################################
##Caminata Aleatoria con y sin drift
set.seed(154)
w = rnorm(200); x = cumsum(w)
wd = w +.2; xd = cumsum(wd)
plot.ts(xd, ylim=c(-5,55), main="Caminata Aletoria", ylab='')
lines(x, col=4); abline(h=0, col=4, lty=2); abline(a=0, b=.2, lty=2)
##### RaÃz unitaria con componentes ARMA
Tlength=200
a0=3
a1=0.5
tiempo=seq(1:Tlength)
xt=a1*tiempo
tendencia=a0+a1*tiempo
drift=a1*rep(1,Tlength)
arimaej=arima.sim(list(order = c(1,0,1), ar = 0.7,ma=0.6), n = Tlength)
plot(arimaej)
caminata=as.ts(cumsum(arimaej))
drift=as.ts(cumsum(arimaej+a0))
#x11()
#par(mfrow=c(2,1))
plot(caminata)
plot(drift) ###Caminata Aleatoria Con Drift
linear=as.ts(cumsum(arimaej+a0+xt))
plot(linear) ###Caminata Aleatoria alrededor de una función cuadrática
auto.arima(arimaej,d=0,D=0,max.p=20,max.q=0,start.p=0, start.q=0,seasonal=FALSE,max.order=20,stationary=TRUE,ic="aic",stepwise=FALSE,allowmean = TRUE)
Series: arimaej
ARIMA(12,0,0) with zero mean
Coefficients:
ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8 ar9 ar10 ar11 ar12
1.2429 -0.7827 0.5695 -0.4683 0.4649 -0.3968 0.3822 -0.3601 0.2715 -0.1908 0.3065 -0.1350
s.e. 0.0700 0.1109 0.1238 0.1293 0.1313 0.1329 0.1322 0.1308 0.1295 0.1262 0.1145 0.0719
sigma^2 = 0.9835: log likelihood = -277.35
AIC=580.71 AICc=582.67 BIC=623.59
ar(arimaej)
Call:
ar(x = arimaej)
Coefficients:
1 2 3 4 5 6 7 8 9
1.2522 -0.7899 0.5314 -0.3955 0.3921 -0.3277 0.3174 -0.2983 0.2001
Order selected 9 sigma^2 estimated as 1.043
fUnitRoots::adfTest(arimaej,lags = 10,type='nc')
Title:
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order: 10
STATISTIC:
Dickey-Fuller: -1.6543
P VALUE:
0.09434
Description:
Wed Aug 28 12:56:56 2024 by user:
aTSA::adf.test(arimaej,nlag = 11)
Augmented Dickey-Fuller Test
alternative: stationary
Type 1: no drift no trend
lag ADF p.value
[1,] 0 -4.02 0.0100
[2,] 1 -5.45 0.0100
[3,] 2 -3.81 0.0100
[4,] 3 -3.77 0.0100
[5,] 4 -3.04 0.0100
[6,] 5 -3.15 0.0100
[7,] 6 -2.72 0.0100
[8,] 7 -2.78 0.0100
[9,] 8 -2.18 0.0298
[10,] 9 -1.96 0.0493
[11,] 10 -1.65 0.0943
Type 2: with drift no trend
lag ADF p.value
[1,] 0 -4.13 0.0100
[2,] 1 -5.68 0.0100
[3,] 2 -3.97 0.0100
[4,] 3 -3.90 0.0100
[5,] 4 -3.14 0.0264
[6,] 5 -3.30 0.0177
[7,] 6 -2.84 0.0584
[8,] 7 -2.92 0.0468
[9,] 8 -2.27 0.2221
[10,] 9 -2.02 0.3185
[11,] 10 -1.70 0.4464
Type 3: with drift and trend
lag ADF p.value
[1,] 0 -4.55 0.0100
[2,] 1 -6.26 0.0100
[3,] 2 -4.50 0.0100
[4,] 3 -4.35 0.0100
[5,] 4 -3.58 0.0362
[6,] 5 -3.87 0.0166
[7,] 6 -3.37 0.0605
[8,] 7 -3.52 0.0419
[9,] 8 -2.81 0.2342
[10,] 9 -2.56 0.3417
[11,] 10 -2.28 0.4566
----
Note: in fact, p.value = 0.01 means p.value <= 0.01
####No hay presencia de RaÃz Unitaria
tseries::adf.test(drift,k = 10) ####Note que esta función no trabaja bien, favor no usar
Augmented Dickey-Fuller Test
data: drift
Dickey-Fuller = -1.9002, Lag order = 10, p-value = 0.6179
alternative hypothesis: stationary
fUnitRoots::adfTest(drift,lags = 10,type='c')
Warning in fUnitRoots::adfTest(drift, lags = 10, type = "c") :
p-value greater than printed p-value
Title:
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order: 10
STATISTIC:
Dickey-Fuller: 1.1031
P VALUE:
0.99
Description:
Wed Aug 28 12:56:57 2024 by user:
aTSA::adf.test(drift,nlag = 11)
Augmented Dickey-Fuller Test
alternative: stationary
Type 1: no drift no trend
lag ADF p.value
[1,] 0 15.74 0.990
[2,] 1 3.13 0.990
[3,] 2 4.20 0.990
[4,] 3 2.89 0.990
[5,] 4 2.94 0.990
[6,] 5 2.57 0.990
[7,] 6 2.64 0.990
[8,] 7 2.38 0.990
[9,] 8 2.41 0.990
[10,] 9 1.89 0.985
[11,] 10 1.75 0.979
Type 2: with drift no trend
lag ADF p.value
[1,] 0 4.76 0.99
[2,] 1 1.58 0.99
[3,] 2 2.13 0.99
[4,] 3 1.45 0.99
[5,] 4 1.53 0.99
[6,] 5 1.48 0.99
[7,] 6 1.52 0.99
[8,] 7 1.44 0.99
[9,] 8 1.46 0.99
[10,] 9 1.16 0.99
[11,] 10 1.10 0.99
Type 3: with drift and trend
lag ADF p.value
[1,] 0 -0.525 0.980
[2,] 1 -1.375 0.837
[3,] 2 -1.080 0.923
[4,] 3 -1.174 0.908
[5,] 4 -1.139 0.914
[6,] 5 -1.439 0.811
[7,] 6 -1.327 0.858
[8,] 7 -1.510 0.781
[9,] 8 -1.420 0.819
[10,] 9 -1.699 0.701
[11,] 10 -1.900 0.616
----
Note: in fact, p.value = 0.01 means p.value <= 0.01
###Otro Ejemplo
estacionario=arima.sim(list(order = c(1,0,1), ar = 0.7,ma=0.6), n = Tlength)
trend.estacionario=tendencia+estacionario
#x11()
plot(trend.estacionario)
fUnitRoots::adfTest(trend.estacionario)
Warning in fUnitRoots::adfTest(trend.estacionario) :
p-value greater than printed p-value
Title:
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order: 1
STATISTIC:
Dickey-Fuller: 2.8525
P VALUE:
0.99
Description:
Wed Aug 28 12:56:57 2024 by user:
fUnitRoots::adfTest(trend.estacionario,lags=1,type='ct') ###Cambie 11 rezago hasta
Warning in fUnitRoots::adfTest(trend.estacionario, lags = 1, type = "ct") :
p-value smaller than printed p-value
Title:
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order: 1
STATISTIC:
Dickey-Fuller: -6.3342
P VALUE:
0.01
Description:
Wed Aug 28 12:56:57 2024 by user:
aTSA::adf.test(trend.estacionario,nlag = 2)
Augmented Dickey-Fuller Test
alternative: stationary
Type 1: no drift no trend
lag ADF p.value
[1,] 0 4.66 0.99
[2,] 1 2.85 0.99
Type 2: with drift no trend
lag ADF p.value
[1,] 0 0.0188 0.957
[2,] 1 -0.2507 0.924
Type 3: with drift and trend
lag ADF p.value
[1,] 0 -4.00 0.0104
[2,] 1 -6.33 0.0100
----
Note: in fact, p.value = 0.01 means p.value <= 0.01
Tlength=200
arimaej_raiz_unit=arima.sim(list(order = c(1,1,1), ar = 0.7,ma=0.6), n = Tlength)
plot(arimaej_raiz_unit)
acf(arimaej_raiz_unit)
pacf(arimaej_raiz_unit)
stats::ar(arimaej_raiz_unit)### Permite Seleccionar el lag
Call:
stats::ar(x = arimaej_raiz_unit)
Coefficients:
1 2 3
1.3077 -0.1686 -0.1512
Order selected 3 sigma^2 estimated as 7.736
fUnitRoots::adfTest(arimaej_raiz_unit,lags =1 ,type='nc')
Title:
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order: 1
STATISTIC:
Dickey-Fuller: -1.2511
P VALUE:
0.2165
Description:
Wed Aug 28 12:57:04 2024 by user:
fUnitRoots::adfTest(arimaej_raiz_unit,type='nc')
Title:
Augmented Dickey-Fuller Test
Test Results:
PARAMETER:
Lag Order: 1
STATISTIC:
Dickey-Fuller: -1.2511
P VALUE:
0.2165
Description:
Wed Aug 28 12:57:04 2024 by user:
####No hay presencia de RaÃz Unitaria
aTSA::adf.test(arimaej_raiz_unit,nlag=12) ####Note que hay
Augmented Dickey-Fuller Test
alternative: stationary
Type 1: no drift no trend
lag ADF p.value
[1,] 0 -0.197 0.587
[2,] 1 -1.251 0.231
[3,] 2 -0.789 0.396
[4,] 3 -1.074 0.294
[5,] 4 -0.841 0.378
[6,] 5 -0.912 0.352
[7,] 6 -0.841 0.378
[8,] 7 -0.975 0.330
[9,] 8 -0.961 0.335
[10,] 9 -0.799 0.393
[11,] 10 -0.796 0.394
[12,] 11 -0.846 0.376
Type 2: with drift no trend
lag ADF p.value
[1,] 0 -1.16 0.641
[2,] 1 -2.11 0.284
[3,] 2 -1.60 0.484
[4,] 3 -1.99 0.330
[5,] 4 -1.63 0.474
[6,] 5 -1.69 0.450
[7,] 6 -1.65 0.464
[8,] 7 -1.75 0.425
[9,] 8 -1.78 0.416
[10,] 9 -1.71 0.441
[11,] 10 -1.66 0.459
[12,] 11 -1.76 0.423
Type 3: with drift and trend
lag ADF p.value
[1,] 0 0.203 0.990
[2,] 1 -2.180 0.499
[3,] 2 -1.331 0.856
[4,] 3 -1.908 0.613
[5,] 4 -1.455 0.804
[6,] 5 -1.568 0.756
[7,] 6 -1.469 0.798
[8,] 7 -1.671 0.713
[9,] 8 -1.664 0.716
[10,] 9 -1.459 0.802
[11,] 10 -1.447 0.807
[12,] 11 -1.539 0.769
----
Note: in fact, p.value = 0.01 means p.value <= 0.01
######
prueba_df=ur.df(arimaej_raiz_unit,type="none",lags=5)
summary(prueba_df)
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-2.41389 -0.72024 0.00503 0.62760 2.77070
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 -0.001714 0.001879 -0.912 0.3628
z.diff.lag1 1.298787 0.072602 17.889 < 2e-16 ***
z.diff.lag2 -0.802310 0.117554 -6.825 1.16e-10 ***
z.diff.lag3 0.533406 0.125282 4.258 3.25e-05 ***
z.diff.lag4 -0.272105 0.117722 -2.311 0.0219 *
z.diff.lag5 0.065239 0.073395 0.889 0.3752
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.004 on 189 degrees of freedom
Multiple R-squared: 0.7757, Adjusted R-squared: 0.7686
F-statistic: 108.9 on 6 and 189 DF, p-value: < 2.2e-16
Value of test-statistic is: -0.9123
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
Vamos a Analizar la Serie de Pasajeros.Iniciamos con las gráficas y transformación de Box-Cox
[1] 4.102259e-05
[1] 0.2
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1949 2.548535 2.561426 2.588461 2.582990 2.567558 2.593773 2.615143 2.615143 2.595510 2.563492 2.529884 2.561426
1950 2.555089 2.577352 2.603953 2.593773 2.575434 2.616686 2.646282 2.646282 2.629992 2.590249 2.552929 2.602296
1951 2.610433 2.618215 2.656336 2.636968 2.648852 2.656336 2.680161 2.680161 2.663501 2.635595 2.612017 2.641022
1952 2.647572 2.658759 2.673698 2.659957 2.662328 2.699069 2.709945 2.720109 2.690390 2.671486 2.648852 2.674793
1953 2.676962 2.676962 2.715111 2.714262 2.709067 2.720927 2.737151 2.742897 2.715955 2.692360 2.658759 2.682259
1954 2.685357 2.668111 2.714262 2.707296 2.713408 2.737151 2.762643 2.756996 2.733443 2.709067 2.684331 2.709067
1955 2.720109 2.712549 2.739332 2.740768 2.741481 2.770427 2.796406 2.787934 2.768668 2.744300 2.715955 2.747066
1956 2.751119 2.746379 2.771588 2.769257 2.772164 2.801154 2.818211 2.814887 2.791987 2.765084 2.742191 2.765084
1957 2.770427 2.762027 2.792485 2.788447 2.791987 2.821853 2.837961 2.838662 2.814466 2.787934 2.764478 2.782161
1958 2.784288 2.772164 2.795437 2.788447 2.795922 2.826939 2.846793 2.851301 2.814466 2.793969 2.767483 2.782696
1959 2.794460 2.785340 2.815307 2.811044 2.821052 2.840400 2.864196 2.867286 2.837255 2.815726 2.795437 2.814887
1960 2.819842 2.808860 2.820650 2.836545 2.840400 2.860440 2.883580 2.879651 2.852247 2.836545 2.808419 2.825783
Ahora avanzamos en el setido de verificar si la serie muestra la presencia de una o varias raÃces unitarias
acf(logAirP,lag.max = 60)
pacf(logAirP)
ar(logAirP) ##Selecciona un modelo AR usando el crierio de Akaike a la serie Aaa
Call:
ar(x = logAirP)
Coefficients:
1 2 3 4 5 6 7 8 9 10 11 12 13
1.0013 -0.0959 0.0329 -0.0252 0.0231 0.0005 -0.0157 -0.0762 0.1071 -0.0236 0.0675 0.4536 -0.4854
Order selected 13 sigma^2 estimated as 0.01336
aTSA::adf.test(logAirP,nlag = 14)
Augmented Dickey-Fuller Test
alternative: stationary
Type 1: no drift no trend
lag ADF p.value
[1,] 0 0.913 0.902
[2,] 1 0.674 0.836
[3,] 2 0.796 0.871
[4,] 3 0.936 0.905
[5,] 4 1.510 0.966
[6,] 5 1.382 0.957
[7,] 6 1.451 0.962
[8,] 7 1.847 0.983
[9,] 8 3.491 0.990
[10,] 9 4.174 0.990
[11,] 10 8.176 0.990
[12,] 11 11.250 0.990
[13,] 12 3.787 0.990
[14,] 13 2.483 0.990
Type 2: with drift no trend
lag ADF p.value
[1,] 0 -1.816 0.400
[2,] 1 -2.018 0.321
[3,] 2 -1.650 0.465
[4,] 3 -1.609 0.481
[5,] 4 -1.288 0.595
[6,] 5 -1.108 0.659
[7,] 6 -0.923 0.723
[8,] 7 -0.807 0.764
[9,] 8 -0.720 0.795
[10,] 9 -0.859 0.746
[11,] 10 -1.578 0.493
[12,] 11 -2.499 0.134
[13,] 12 -1.952 0.347
[14,] 13 -1.717 0.439
Type 3: with drift and trend
lag ADF p.value
[1,] 0 -4.85 0.0100
[2,] 1 -7.00 0.0100
[3,] 2 -6.71 0.0100
[4,] 3 -7.13 0.0100
[5,] 4 -5.66 0.0100
[6,] 5 -6.42 0.0100
[7,] 6 -6.88 0.0100
[8,] 7 -6.28 0.0100
[9,] 8 -3.62 0.0341
[10,] 9 -2.98 0.1692
[11,] 10 -1.32 0.8577
[12,] 11 -0.44 0.9833
[13,] 12 -1.53 0.7696
[14,] 13 -2.15 0.5109
----
Note: in fact, p.value = 0.01 means p.value <= 0.01
summary(ur.df(logAirP,type="none",lags = 13))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-0.090906 -0.022763 -0.002391 0.022181 0.107123
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 0.003636 0.001464 2.483 0.014448 *
z.diff.lag1 -0.379495 0.088572 -4.285 3.8e-05 ***
z.diff.lag2 -0.218029 0.071925 -3.031 0.003003 **
z.diff.lag3 -0.137513 0.074265 -1.852 0.066618 .
z.diff.lag4 -0.227017 0.073095 -3.106 0.002386 **
z.diff.lag5 -0.098892 0.075307 -1.313 0.191712
z.diff.lag6 -0.193729 0.071366 -2.715 0.007649 **
z.diff.lag7 -0.153409 0.072344 -2.121 0.036090 *
z.diff.lag8 -0.264234 0.072215 -3.659 0.000382 ***
z.diff.lag9 -0.100257 0.075728 -1.324 0.188138
z.diff.lag10 -0.204136 0.072692 -2.808 0.005847 **
z.diff.lag11 -0.088297 0.073995 -1.193 0.235190
z.diff.lag12 0.694627 0.072012 9.646 < 2e-16 ***
z.diff.lag13 0.311279 0.088880 3.502 0.000656 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.04121 on 116 degrees of freedom
Multiple R-squared: 0.8701, Adjusted R-squared: 0.8544
F-statistic: 55.5 on 14 and 116 DF, p-value: < 2.2e-16
Value of test-statistic is: 2.4833
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
summary(ur.df(logAirP,type="drift",lags = 13))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression drift
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-0.097948 -0.021365 -0.003174 0.022461 0.101422
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.120185 0.056858 2.114 0.036696 *
z.lag.1 -0.016718 0.009737 -1.717 0.088668 .
z.diff.lag1 -0.401952 0.087921 -4.572 1.23e-05 ***
z.diff.lag2 -0.253986 0.072886 -3.485 0.000698 ***
z.diff.lag3 -0.180754 0.075984 -2.379 0.019015 *
z.diff.lag4 -0.266718 0.074435 -3.583 0.000499 ***
z.diff.lag5 -0.143868 0.077196 -1.864 0.064919 .
z.diff.lag6 -0.234797 0.072957 -3.218 0.001676 **
z.diff.lag7 -0.199130 0.074495 -2.673 0.008609 **
z.diff.lag8 -0.310371 0.074432 -4.170 5.94e-05 ***
z.diff.lag9 -0.155320 0.079037 -1.965 0.051808 .
z.diff.lag10 -0.255366 0.075618 -3.377 0.001000 **
z.diff.lag11 -0.143773 0.077492 -1.855 0.066112 .
z.diff.lag12 0.641419 0.075291 8.519 7.14e-14 ***
z.diff.lag13 0.290254 0.088144 3.293 0.001318 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.04061 on 115 degrees of freedom
Multiple R-squared: 0.874, Adjusted R-squared: 0.8587
F-statistic: 56.98 on 14 and 115 DF, p-value: < 2.2e-16
Value of test-statistic is: -1.717 5.4095
Critical values for test statistics:
1pct 5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1 6.52 4.63 3.81
######SE debe chequear si hay que diferenciar de nuevo!!!
dlogAirPass=diff(logAirP)
plot(dlogAirPass)
#####Transformación requrida para los datos(transformación Box-Cox y diferencia ordinaria)
pacf(dlogAirPass)
acf(dlogAirPass,lag.max = 48)
ar(dlogAirPass)
Call:
ar(x = dlogAirPass)
Coefficients:
1 2 3 4 5 6 7 8 9 10 11 12 13
-0.0993 -0.1226 -0.2824 -0.2875 -0.1184 -0.2406 -0.1829 -0.2971 -0.1040 -0.2627 -0.1259 0.5734 0.0179
14 15
-0.1667 0.1200
Order selected 15 sigma^2 estimated as 0.002988
aTSA::adf.test(dlogAirPass,nlag = 15)
Augmented Dickey-Fuller Test
alternative: stationary
Type 1: no drift no trend
lag ADF p.value
[1,] 0 -9.606 0.0100
[2,] 1 -8.816 0.0100
[3,] 2 -7.634 0.0100
[4,] 3 -8.745 0.0100
[5,] 4 -6.791 0.0100
[6,] 5 -6.245 0.0100
[7,] 6 -6.571 0.0100
[8,] 7 -9.498 0.0100
[9,] 8 -8.028 0.0100
[10,] 9 -10.214 0.0100
[11,] 10 -7.327 0.0100
[12,] 11 -1.868 0.0622
[13,] 12 -1.223 0.2402
[14,] 13 -1.283 0.2186
[15,] 14 -0.992 0.3231
Type 2: with drift no trend
lag ADF p.value
[1,] 0 -9.63 0.0100
[2,] 1 -8.86 0.0100
[3,] 2 -7.71 0.0100
[4,] 3 -8.94 0.0100
[5,] 4 -6.98 0.0100
[6,] 5 -6.46 0.0100
[7,] 6 -6.91 0.0100
[8,] 7 -10.55 0.0100
[9,] 8 -9.57 0.0100
[10,] 9 -15.32 0.0100
[11,] 10 -16.05 0.0100
[12,] 11 -4.60 0.0100
[13,] 12 -3.05 0.0352
[14,] 13 -3.22 0.0222
[15,] 14 -2.72 0.0789
Type 3: with drift and trend
lag ADF p.value
[1,] 0 -9.60 0.0100
[2,] 1 -8.83 0.0100
[3,] 2 -7.69 0.0100
[4,] 3 -8.92 0.0100
[5,] 4 -6.95 0.0100
[6,] 5 -6.43 0.0100
[7,] 6 -6.88 0.0100
[8,] 7 -10.51 0.0100
[9,] 8 -9.55 0.0100
[10,] 9 -15.44 0.0100
[11,] 10 -16.57 0.0100
[12,] 11 -4.94 0.0100
[13,] 12 -3.37 0.0628
[14,] 13 -3.56 0.0392
[15,] 14 -3.12 0.1090
----
Note: in fact, p.value = 0.01 means p.value <= 0.01
summary(ur.df(dlogAirPass,type="none",lags = 14))
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################
Test regression none
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
Residuals:
Min 1Q Median 3Q Max
-0.081539 -0.022167 0.001427 0.028987 0.098928
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1 -0.33366 0.33626 -0.992 0.323188
z.diff.lag1 -0.94573 0.33797 -2.798 0.006042 **
z.diff.lag2 -1.06221 0.34615 -3.069 0.002692 **
z.diff.lag3 -1.22917 0.34971 -3.515 0.000634 ***
z.diff.lag4 -1.32684 0.32990 -4.022 0.000105 ***
z.diff.lag5 -1.26830 0.30882 -4.107 7.61e-05 ***
z.diff.lag6 -1.34019 0.28801 -4.653 8.95e-06 ***
z.diff.lag7 -1.33665 0.26537 -5.037 1.81e-06 ***
z.diff.lag8 -1.46466 0.24807 -5.904 3.79e-08 ***
z.diff.lag9 -1.41286 0.23488 -6.015 2.26e-08 ***
z.diff.lag10 -1.49206 0.22124 -6.744 6.87e-10 ***
z.diff.lag11 -1.42586 0.20715 -6.883 3.46e-10 ***
z.diff.lag12 -0.60007 0.19564 -3.067 0.002704 **
z.diff.lag13 -0.22366 0.15295 -1.462 0.146431
z.diff.lag14 -0.20779 0.09268 -2.242 0.026919 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.04183 on 113 degrees of freedom
Multiple R-squared: 0.9165, Adjusted R-squared: 0.9054
F-statistic: 82.71 on 15 and 113 DF, p-value: < 2.2e-16
Value of test-statistic is: -0.9923
Critical values for test statistics:
1pct 5pct 10pct
tau1 -2.58 -1.95 -1.62
monthplot(dlogAirPass)
spectrum(dlogAirPass)
Estimated spectral density
series: dlogAirPass
method: Raw Periodogram
nseasons: 12
frequency range: (-6,6]
sort method for the table: decreasing magnitudes
freq spec % Total Cum. %
[1,] 2.000000 1.73274552 0.28973097 0.2897310
[2,] 1.000000 1.61152998 0.26946262 0.5591936
[3,] 4.000000 0.92289692 0.15431684 0.7135104
[4,] 5.000000 0.48378566 0.08089341 0.7944038
[5,] 3.000000 0.46145560 0.07715962 0.8715635
[6,] 4.166667 0.08107252 0.01355607 0.8851195
####Identificación de los Órdenes Autoregresivos
acf(dlogAirPass,lag.max = 60,ci.type='ma') ###Se requiere un MA de orden muy grande
pacf(dlogAirPass,lag.max = 60) ####Puede ser un autoregresivo de orden 12
###Arima Automático
modelo.automatico1=auto.arima(dlogAirPass,d=0,D=0,max.p=12,max.q=48,start.p=0, start.q=0,seasonal=FALSE,max.order=12,stationary=TRUE,ic="aicc",stepwise=FALSE,allowmean = TRUE)
#####Ajuste del Modelo
####Note que entramos la serie original
library(TSA)
Attaching package: ‘TSA’
The following objects are masked from ‘package:PerformanceAnalytics’:
kurtosis, skewness
The following object is masked from ‘package:readr’:
spec
The following object is masked from ‘package:sarima’:
periodogram
The following objects are masked from ‘package:stats’:
acf, arima
The following object is masked from ‘package:utils’:
tar
AjusteArima93=forecast::Arima(AirPassengers,order = c(9,1,3),lambda = 0,include.drift = TRUE)
summary(AjusteArima93)
Series: AirPassengers
ARIMA(9,1,3) with drift
Box Cox transformation: lambda= 0
Coefficients:
ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8 ar9 ma1 ma2 ma3
-0.4424 -0.2037 0.0860 -0.6827 -0.271 -0.1659 -0.1491 -0.6805 -0.2141 0.4821 -0.1477 -0.7510
s.e. 0.1105 0.0758 0.0759 0.0713 0.100 0.0711 0.0750 0.0715 0.1037 0.0821 0.0677 0.0689
drift
0.0101
s.e. 0.0009
sigma^2 = 0.004642: log likelihood = 183.78
AIC=-339.55 AICc=-336.27 BIC=-298.07
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.2192766 20.40352 15.27227 -0.2117784 5.395308 0.476807 0.012135
coeftest(AjusteArima93)
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 -0.44238067 0.11049430 -4.0037 6.237e-05 ***
ar2 -0.20368876 0.07583340 -2.6860 0.007231 **
ar3 0.08599219 0.07591330 1.1328 0.257311
ar4 -0.68273511 0.07128722 -9.5772 < 2.2e-16 ***
ar5 -0.27099739 0.10000994 -2.7097 0.006734 **
ar6 -0.16591360 0.07107162 -2.3345 0.019572 *
ar7 -0.14906334 0.07498938 -1.9878 0.046835 *
ar8 -0.68051426 0.07152435 -9.5144 < 2.2e-16 ***
ar9 -0.21410059 0.10371315 -2.0644 0.038984 *
ma1 0.48207457 0.08205276 5.8752 4.224e-09 ***
ma2 -0.14767013 0.06771628 -2.1807 0.029204 *
ma3 -0.75098362 0.06886761 -10.9047 < 2.2e-16 ***
drift 0.01009809 0.00090294 11.1836 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#####Refinando el modelo
AjusteArima12=forecast::Arima(AirPassengers,order = c(12,1,0),lambda = 0,include.drift = TRUE)
#,fixed=c(NA,NA,0,NA,NA,NA,NA,NA,NA,NA,0,NA,NA)
summary(AjusteArima12)
Series: AirPassengers
ARIMA(12,1,0) with drift
Box Cox transformation: lambda= 0
Coefficients:
ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8 ar9 ar10 ar11 ar12
-0.2212 -0.2835 -0.2550 -0.2973 -0.2299 -0.2622 -0.257 -0.3274 -0.2194 -0.2934 -0.198 0.6020
s.e. 0.0670 0.0679 0.0673 0.0699 0.0679 0.0657 0.066 0.0679 0.0684 0.0681 0.068 0.0683
drift
0.0098
s.e. 0.0011
sigma^2 = 0.001807: log likelihood = 246.66
AIC=-465.31 AICc=-462.03 BIC=-423.83
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.1137697 10.82763 8.086473 -0.007108351 3.093656 0.2524632 -0.2063749
coeftest(AjusteArima12)
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 -0.2212162 0.0669842 -3.3025 0.0009582 ***
ar2 -0.2834744 0.0679448 -4.1721 3.018e-05 ***
ar3 -0.2549912 0.0672627 -3.7910 0.0001501 ***
ar4 -0.2973403 0.0698594 -4.2563 2.079e-05 ***
ar5 -0.2298510 0.0679217 -3.3841 0.0007142 ***
ar6 -0.2621654 0.0657330 -3.9883 6.654e-05 ***
ar7 -0.2570036 0.0660119 -3.8933 9.889e-05 ***
ar8 -0.3274286 0.0679449 -4.8190 1.443e-06 ***
ar9 -0.2194360 0.0684078 -3.2078 0.0013377 **
ar10 -0.2934215 0.0680902 -4.3093 1.638e-05 ***
ar11 -0.1979913 0.0680155 -2.9110 0.0036031 **
ar12 0.6020203 0.0682917 8.8154 < 2.2e-16 ***
drift 0.0097712 0.0010737 9.1003 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
residuales=AjusteArima12$residuals
plot(residuales)
#
#plot(SDresiduales)
acf(residuales,lag.max = 48)
pacf(residuales)
#acf(SDresiduales)
#pacf(SDresiduales)
#Test de normalidad
jarque.bera.test(residuales)
Jarque Bera Test
data: residuales
X-squared = 7.0597, df = 2, p-value = 0.02931
#Test de autocorrelación
length(residuales)/4
[1] 36
sqrt(length(residuales))
[1] 12
Box.test(residuales, lag =36 , type = "Ljung-Box", fitdf = 12)
Box-Ljung test
data: residuales
X-squared = 64.81, df = 24, p-value = 1.298e-05
monthplot(residuales)
spectrum(residuales,spans = c(3,5))
Estimated spectral density
series: residuales
method: Smoothed Periodogram
nseasons: 12
frequency range: (-6,6]
sort method for the table: decreasing magnitudes
freq spec % Total Cum. %
[1,] 5.583333 0.2454543 0.04105811 0.04105811
[2,] 5.666667 0.2217983 0.03710108 0.07815919
[3,] 5.500000 0.2216797 0.03708124 0.11524043
[4,] 5.416667 0.1780115 0.02977668 0.14501711
[5,] 5.750000 0.1645269 0.02752105 0.17253816
[6,] 2.666667 0.1408333 0.02355775 0.19609590
###Estad?sticas CUSUM
res=residuales
cum=cumsum(res)/sd(res)
N=length(res)
cumq=cumsum(res^2)/sum(res^2)
Af=0.948 ###Cuantil del 95% para la estad?stica cusum
co=0.14013####Valor del cuantil aproximado para cusumsq para n/2
####Para el caso de la serie de pasajeros es aprox (144-12)/2=66
LS=Af*sqrt(N)+2*Af*c(1:length(res))/sqrt(N)
LI=-LS
LQS=co+(1:length(res))/N
LQI=-co+(1:length(res))/N
plot(cum,type="l",ylim=c(min(LI),max(LS)),xlab="t",ylab="",main="CUSUM")
lines(LS,type="S",col="red")
lines(LI,type="S",col="red")
#CUSUMSQ
plot(cumq,type="l",xlab="t",ylab="",main="CUSUMSQ")
lines(LQS,type="S",col="red")
lines(LQI,type="S",col="red")
#####Fase de Pronósticos
pronosticos12=forecast::forecast(AjusteArima12,h=12,level=0.95)
plot(pronosticos12)